Quantum Computing of Poincare Recurrences and Periodic Orbits 
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Quantum algorithms are built enabling to find Poincare recurrence times and periodic orbits of 
classical dynamical systems. It is shown that exponential gain compared to classical algorithms 
can be reached for a restricted class of systems. Quadratic gain can be achieved for a larger set of 
dynamical systems. The simplest cases can be implemented with small number of qubits. 



PACS numbers: 03.67.Lx, 05.45.Ac, 05.45.Tp 

It has been suggested since Feynman [1] that the super- 
position principle of quantum mechanics enables to per- 
form exponentially many computations in parallel (see 
reviews in [2-4]). Thus in principle quantum processors 
using the full power of quantum mechanics may be enor- 
mously faster than classical nowadays computers. This 
possibility has motivated a great deal of attention in the 
scientific community, and many experimental proposals 
are explored to realize such a quantum computer. How- 
ever, it has been surprisingly hard to spot specific prob- 
lems where quantum algorithms may be faster than clas- 
sical ones. One such example is the celebrated Shor's al- 
gorithm [5] which factors large numbers with exponential 
efficiency compare to any known classical algorithm. An- 
other algorithm, due to Grover [6], speeds up the search 
in an unsorted database, although not exponentially. 

In parallel, people have investigated the possibility of 
using a quantum device to simulate physical systems, a 
problem of much practical interest. In particular, several 
works have shown that quantum computers can speed up 
the simulation of quantum mechanical systems [7-9] and 
classical spin systems [10]. It was also suggested in [11] 
that a quantum computer may be efficient at simulat- 
ing classical dynamical systems as well. In the present 
work, the problem of finding Poincare recurrence times 
and periodic orbits of dynamical systems is studied. Such 
quantities give important information about the dynam- 
ics, and are standard tools for studying classical mechan- 
ics. Explicit algorithms will be presented showing that 
such quantities can be evaluated faster on a quantum 
computer than on a classical one. In some cases, expo- 
nential efficiency compared to classical algorithms can be 
reached. In other cases, quadratic gains can be realized. 

For classical bounded conservative systems, Poincare 
[12] showed that some points from an arbitrary small 
phase space domain will eventually come back to this do- 
main. In general, actual values of the recurrence time for 
a given area of phase space, and statistics of these times, 
are not given by the theorem and require specific studies 
to be understood. This time can be extremely long, and 
is very hard to study numerically. Advances in the power 
of modern classical computers allow to get information 



about such quantities, but the problem is still under in- 
vestigation nowadays [13]. In addition to their intrinsic 
interest, the return times give also insight on transport 
properties, diffusion coefficients and correlation functions 
of the system studied [13]. A related feature of dynami- 
cal systems is the presence of periodic orbits, which are 
trajectories coming back exactly to their initial point in 
phase space. They play a special role in understand- 
ing physical properties of these systems, and are present 
both in Hamiltonian and dissipative systems. For ex- 
ample, they enable to compute diffusion coefficients and 
properties of strange attractors [14], and enter the famous 
Gutzwiller trace formula [15] which is a general semiclas- 
sical quantization scheme valid for chaotic systems. 

I will first show how such quantities can be found ef- 
ficiently on a particular example, the Arnold cat map 
[16,17], one of the most famous chaotic dynamical sys- 
tems. It is an automorphism of the torus with equation: 



y = y + x (mod 1) , x = y + 2x (mod 1) 



(1) 



where bars denote the new values of the variables after 
one iteration. This map is area-preserving, and may be 
interpreted as describing the evolution in phase space, 
with x representing position and y momentum. It is 
an Anosov system, with homogeneous exponential diver- 
gence of trajectories and positive Kolmogorov-Sinai en- 
tropy h rs 0.96. The map (1) can be equivalently written 

as action of the 2x2 matrix L 
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Other cat maps can be built by taking any 2x2 matrix 
L with integer entries and determinant 1. 

It is known [18-20] that for the cat map (1) periodic 
orbits are in one-to-one correspondence with points with 
rational coordinates. Still, it is a computationally hard 
problem to find the period of a given point. An usual 
way for studying this problem is to consider the set of 
rational points sharing the same denominator g. These 
points form a g x g lattice in the phase space, which is 
invariant under the action of the map. The action of the 
map on such points can be written on the numerators 
only, namely: y — y + x (mod g) , x — y + 2x (mod g) , 
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or ^ _ j = L j (mod g) , with x, y, x, y integers. 

All points in such a lattice are periodic, but not with 
the same period. A quantity similar to recurrence time 
is the lattice period function a(g). It is the smallest 
integer such that after a(g) iterations all points in the 
lattice have come back to the initial position, i.e. it is 
the recurrence time of the whole lattice. This quantity 
has been the subject of many studies [18-20], since it 
describes the whole set of periodic orbits of the system, 
and also gives insight on the quantization of the cat map 
[18,20]. In particular, the behavior of a{g) controls the 
crgodicity of the eigenfunctions of the quantum cat map 
for Ti = 1/g [21]. It turns out to be a very erratic func- 
tion of g; some results on its statistical properties have 
been derived nonrigorously based on assumptions from 
probabilistic number theory [20]. Numerical checks for 
moderate values of g show that the asymptotic regime is 
very slowly reached [20]. Thus for example computing 
a(g) in small windows around large values of g will en- 
able to check its asymptotic properties. This has interest 
both for number theory and study of dynamical systems. 

A quantum computer of moderate size is able to com- 
pute the function a(g) for exponentially large values of 
g. One does not need to actually simulate the evolu- 
tion of individual points in the lattice, since a(g) is the 
smallest integer t such that L* = I (mod g) , with / the 
identity. Nine registers of n q qubits with n q ~ log 2 g are 
needed, one holding the values of time, the next four the 
matrix entries, and the last ones being workspace. The 
algorithm applies period-finding (which is the basis of Si- 
mon's [22] or Shor's [5] algorithms) to the four registers 
encoding the entries of the matrix. First one prepares the 
initial state, which is A" 1 / 2 ]^^ 1 l^l 1 )! )! )! 1 ) where 
N = 2 n ". It is easily built from the ground state of 
the system by applying n q Hadamard gates and 2 singlc- 
qubit flips. Then one constructs the time evolution to 
end up with the entries of the matrix L* modulo g in 
the last four registers. To this aim, first the entries 
(di,bi,Ci,di) of each L 2% for i — l,...,n q — 1 are pre- 
computed classically. This can be done sequentially, first 
squaring L to get L 2 , then squaring L 2 and so on until 
L 2 9 is obtained. This requires of the order of 0(n q ) 
classical operations. For each t, the binary decompo- 
sition of t is t = £™! oii2 l with a>i = 0,1, so that 
L* = n™! (L 2 ) ai . So the exponentiation of L can be 
done by n q multiplication by L 2 conditioned on the value 
of the qubit i. This matrix multiplication can be done by 
a sequence of number multiplication. For example, the 
sequence of transformations for the first two entries is: 

iV-i/2 E ;y |t)|a)|6)|c)|d)|a)|6) (copy) 

- N- 1 / 2 ^- 1 ^aa^bd^c^ab^ba) 

-» AT-^Eto 1 + bc i)\ ab i + bdi)\c)\d)\abi)\bci). 

This needs a controlled multiplier modulo g, which can 



be done following the procedure in [23] using two ex- 
tra registers as workspace. To erase the unwanted reg- 
isters one uses a standard trick: one builds the se- 
quence of gates corresponding to F : |u)|u)|0)|0) — > 
\u)\v)\dibiU — Cibiv)\ciiCiV — biau). Applying F~ x to 
the state above sets the work registers to zero. One 
then does again the same operation for the next two en- 
tries of the matrix. This needs 0(n 2 ) operations, and 
should be done for all qubits i in the first register to 
build L*, so 0(n q ) gates are used in total. One ends 
up with the state A" 1 / 2 ^^q 1 \t)\A t )\B t )\C t )\D t ) where 
(A t , B t , C t ,D t ) are entries of the matrix L* modulo g. 

Then as in Shor's algorithm one performs a Quantum 
Fourier Transform (QFT) of the first register. This quan- 
tum analog of the classical Fourier transform needs 0(n 2 ) 
elementary gates to be performed. Since the last four 
registers describe a periodic function of t, peaks in the 
Fourier transform can be measured and yield the period 
ct(g). Indeed, if a(g) divides A, then peaks are precisely 
at values of multiples of N/a(g), while otherwise broader 
peaks will appear but significant probability will be con- 
centrated on the best rational approximants of multiples 
of N/a(g), and a continued fractions algorithm can yield 
the period in 0(n q ) classical operations [2-4]. 

It is worth trying to compare this algorithm to com- 
pute a(g) with others. Classically, computing all iterates 
L k up to the period is clearly exponential since a(g) is of 
the same order of magnitude as g ((a(g)/g) — > slower 
than any inverse power of g [20]). Another method to 
compute a(g) [19] is based on the theory of ideals of 
quadratic fields. For a prime p, a(p) = p [^/^ , where 
() is the Kronecker symbol (= ±1), d is the discriminant 
of the field containing the eigenvalues of L, and m is an 
integer divisor of p— (d/p) to be determined by trial and 
errors. For composite g, a(g) can be determined from the 
a(p), p being the prime factors of g. This method is also 
exponential classically, since it requires the factorization 
of g in prime factors. On a quantum computer, one may 
try to use Shor's algorithm, first factoring g, then fac- 
toring p — (d/p) for all prime factors of g and trying all 
possible integer divisors. Most integers g have ~ log log g 
prime factors and for most p there are ~ (logp) log2 divi- 
sors of p— (d/p) [24]. Thus this algorithm is also poly- 
nomial for most numbers. Still, it is known [24] that the 
number of divisors of an integer n can be quite large, 
of the order of 2 log "/ log log ™ in worst cases, making this 
method not polynomial for some g. The algorithm here is 
simpler, giving a(g) in one run, and always polynomial. 

Another possibility to explore periodic orbits and 
Poincare recurrences which can be generalized to other 
systems beyond (1), is to start from an initial point and 
look for periodicities of its iterates. In this case, it is ap- 
propriate to use a discretized map. As discussed in [11], 
the best possible discretization uses symplectic maps [25] . 
These maps can be exactly iterated and are themselves 
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Hamiltonian if the original system is. In dimension 2, 
they can be written as the action of a unitary operator L 
on the N 2 points yj) with Xi — i/N, i = 0, ...,N — 1 
and yj = j/N, j = 0, ...,N - 1, with N — 2 n " . For ex- 
ample, for the cat map, a point (x^ j/j) of the discretized 
phase space density is mapped to (x~i,y~j) through the ac- 
tion of the 2x2 matrix L defined above. The algorithm 
requires three registers, two of size n q specifying points in 
phase space and one of size p holding the values of time, 
plus additional qubits as workspace, with usually p w n q . 

One starts from the initial state 2~ p / 2 X^o* \t)\ x o)\Do)i 
easily built from the action of p Hadamard gates and 
2n q single-qubit flips. Then one constructs a sequence 

of gates giving 2-p/ 2 ^- 1 |t> j X* (a;o)> (2/o)> ■ Then a 
QFT of the first register will yield peaks, from which the 
period of the point can be found. If L and each L 2 
for k = l,...,p— 1 can be computed and implemented 
in (0(B(n q )), where B is a polynomial, this will be ex- 
ponentially faster than classically. For the cat map, all 
can be done efficiently by a method similar to the one 
described above, first precomputing classically the L 2 , 
then decomposing each t in binary representation, and 
then using a controlled multiplier modulo N to compute 
(L t (xo),L t (yo)). Again, this uses 0(n q ) quantum gates. 

Other maps can be explored in the same way, pro- 
vided their classical evolution operator L enables fast 
(polynomial) classical computation of (L 2 ). In such 
cases, the algorithm above can reach exponential gain 
over the classical computation. If L is efficiently implc- 
mentable, but exponentially large iterates cannot be ef- 
ficiently computed (see examples below) the algorithm 
above will not work. Nevertheless, in this case it is 
possible to obtain a gain for the computation of peri- 
odic orbits and recurrence properties of the whole system 
for a given time (as opposed to a chosen initial point). 
The basic idea is to compute sequentially the iterates of 
a distribution and then use Grover iterations to search 
for specific trajectories. The dynamical system is dis- 
cretized through a symplectic map L, on a lattice of 
size N x N where N = 2 n " , and it is assumed that L 
can be implemented on these N 2 points in 0(B(n q )) 
operations where B is a polynomial. Then a subdo- 
main A is selected. The simplest case is a square of size 
P x P with P = 2 P and p < n q . The initial state is 
=o \ x i)\Vi)- Then one realizes the 
transformation |^ ) -» 2~ p Ya=o 1 l L ( x 0> \ L (Uj)) 

a-'E^o'Ejlo 1 |i*(^i)>|i*(l/i)>- The value of the 
first n q — p qubits of each register is enough to know 
if the trajectory is back to the domain A (since A can 
be moved efficiently to a corner). So after t iterations 
this value is checked for all trajectories and a a z con- 
trolled by the value of these 2n q — 2p qubits gives a 
minus sign to the values of \L t (xi))\L t (yj)) which end 
a trajectory which returns to A. This can be done ef- 



ficiently using 0((2n q — 2p) 2 ) elementary gates and no 
work qubit [4]. Then one inverts all the gates but the 

last operation, ending with 2~ p J2i=o Sj=o e \ x i)\Vj}> 
where e = ±1 depending on the fact that the trajec- 
tory went back to A after these t iterations. This whole 
procedure is then used as an oracle for iterations of the 
Grover algorithm [6]. After P/^/M iterations, where M 
is the number of solutions, the amplitude of the wave- 
function is concentrated on initial points which return to 
A after t iterations, and a measure of the registers will 
give one of them. The whole algorithm necessitates 2n q 
qubits plus the workspaces. Only 0(tP/vM) operations 
are needed to get the returns, as opposed to 0(tP 2 /M) 
for the classical simulation (up to logarithmic factors). 
In 0(tP j \J~M) one initial point of a trajectory which re- 
turns to A is obtained. To get other trajectories which 
return, or for different return times, one needs to restart 
the algorithm. This process has obviously interest mostly 
for t <C P. The number of returning trajectories M can 
be evaluated through the quantum counting algorithm 
(phase estimation combined with the Grover iterations) 
in 0(tP) operations [26]. Since generally the number M 
is unknown, one can estimate it by quantum counting 
prior to the search, but it is not necessary [26]. 

The algorithm can be modified to get a differ- 
ent result, which is the number and coordinates of 
the periodic orbits of the system. In this case 
one starts from all the N x N points of the lat- 
tice with N = 2 n ", with initial state |^o) = 

2_n, Z)i=o _ Sj=cT \ x i)\Vj)- Tnen onc performs: 

i^oh 2-»«Ero" i E;io _i ki>iw>ii(^)>ii(w)> 

a-^ECo^EjCo" 1 IxMvMVixiMLtfa)). After t it- 
erations the value of the iterate is compared to the initial 
value (for example by adding the registers bitwise modulo 
2), and a minus sign is given (by using a a z controlled by 
the value of these qubits [4]) if it is the same. Then one 
invert all the gates but the last operation, ending with 

2_ " , Ei=o _1 Ej=o _le l x «)l%>' where e = ±1 depending 
on the fact that the trajectory is periodic of period t or 
not. This whole procedure is then used as an oracle for it- 
erations of the Grover algorithm [6] and all periodic orbits 
can be found by iterating the process. After 0(tN/\/M) 
operations, the amplitude of the wavefunction is concen- 
trated on periodic points of period t (M is the number 
of such points). In contrast, the classical search needs 
0(tN 2 /M) operations. As above M can be evaluated 
through quantum counting in 0(tN) operations. 

A much studied class of classical twist maps corre- 
sponds to the form n = n — KV' (9)(mod 2ttL); 6 = 
9 + n(mod 2ir) where (n, 9) is the pair of conjugated 
momentum (action) and angle variables, and the bars de- 
note the resulting variables after one iteration of the map. 
The discretized map on a N x N lattice is simply Y — 
Y +[NKV'(27rX/N)/(2Tr)]{modN); X = A + YXmodAO 
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where [...] is the integer part and X, Y are integers. The 
case V{9) — cos 9 corresponds to the Chirikov stan- 
dard map which has been a cornerstone model in the 
study of chaos [17]. The discretized map on a 2 nq x 2 Uq 
lattice can be implemented in O(n^) gates on a quan- 
tum computer [8] and therefore a quadratic gain can be 
achieved for return times and periodic orbits (although 
other quantities may be obtained with exponential gain 
[11]). A simpler example is the sawtooth map which 
corresponds to V{9) = —9 2 /2. The discretized map- 
ping is Y = Y + [NK(2ttX/N - 7r)/(27r)](modA); X = 
X + y(modiV). This model has been intensively studied 
in the field of classical chaos [27] . Depending on the val- 
ues of K, the system can be stable or chaotic, and can dis- 
play complex structures in phase space with chaotic and 
integrable parts. Interesting phenomena such as anoma- 
lous diffusion, self-similar structures, etc... are present. 
For integer K, exponentially large iterates can be com- 
puted efficiently and the gain is exponential. For non in- 
teger K, the iterates cannot be computed efficiently but 
the map itself can be implemented in 0(n q ) gates on a 
2"« x 2 Uq lattice, and a quadratic gain can be achieved for 
return times and periodic orbits. For example, the case 
K = ±1/2 for return times requires only three registers 
(two for X and Y, one as workspace) and the algorithm 
can be realized with as little as 8 qubits and 40 gates per 
Grover iteration (domain 4x4ina8x8 lattice). 

A particular case arises for systems with symmetries. 
For example, in the standard and sawtooth maps, the 
classical evolution operator L is the product of two invo- 
lutions L = with l\ = /| = 1. Such involutions have 
lines of fixed points, and it is easy to show that trajecto- 
ries crossing twice one of these lines are periodic. They 
correspond to periodic orbits which are invariant through 
one of the symmetries. These orbits can be found by it- 
crating only points from one of these lines, i.e. N points 
instead of N 2 . The latter algorithm finds these orbits in 
y/~N operations, keeping the quadratic speed-up. 

The periodic orbits found by the various algorithms 
above are exact periodic orbits of the discretized sys- 
tems, which is a Hamiltonian system in its own right. As 
concerns the original continuous system, for hyperbolic 
systems the shadowing theorem [28] ensures that an ex- 
act trajectory will remain close to the dynamics of each 
discretized point for arbitrary times. This exact trajec- 
tory will return close to its starting point after the time 
computed by the algorithm, meaning that it is a Poincare 
recurrence time of the system. Moreover, such approxi- 
mately periodic trajectory can be used as a starting point 
for a Newton method converging to a true periodic orbit. 

In conclusion, it has been shown that on a quantum 
computer one can obtain Poincare recurrence times and 
periodic orbits of certain classical dynamical systems ex- 
ponentially faster that on a classical computer. For a 
larger class of systems, a quadratic gain can be achieved. 
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